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Abstract. Hard X-ray spectra in solar flares provide knowledge of the electron 
spectrum that results from acceleration and propagation in the solar atmosphere. 
However, the inference of the electron spectra from solar X-ray spectra is an ill-posed 
inverse problem. Here we develop and apply an enhanced regularization algorithm 
for this process making use of physical constraints on the form of the electron 
spectrum. The algorithm incorporates various features not heretofore employed in 
the solar flare context : Generalized Singular Value Decomposition (GSVD) to deal 
with different orders of constraints; rectangular form of the cross-section matrix to 
extend the solution energy range; regularization with various forms of the smoothing 
operator; and "preconditioning" of the problem. We show by simulations that this 
technique yields electron spectra with considerably more information and higher 
quality than previous algorithms. 
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1. Introduction 



In order to address fundamental questions on electron propagation and 
acceleration in solar flares, it is necessary to infer as much as quan- 
titative information as possible on the electron spectrum in the solar 
plasma. A longstanding method of doing this involves the analysis of 
the emitted hard X-ray (HXR) bremsstrahlung spectrum, in particular 
the inversion of the integral equation (Brown, 1971) relating the two 
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spectra. This task is particularly challenging since even the most accu- 
rate photon spectra are contaminated by noise, which is dramatically 
amplified in any unconstrained attempt to extract the electron flux 
spectrum. Traditional approaches to the determination of the electron 
spectrum sidestep this problem by assuming simple (e.g. isothermal + 
power-law) forms and adjusting their parameters to achieve the best fit 
to the HXR data (e.g., Holman et al. 2003). However, such algorithms, 
by their very nature, cannot detect features in the electron spectrum 
that, although real, were not included into the prescribed empirical 
form. Indeed, the analysis of high resolution RHESSI (Ramaty High 
Energy Solar Spectroscopic Imager) (Lin et al., 2002) spectra shows 
substantial deviations from simple models (Kontar et. al., 2003). While 
this was adequate for earlier low resolution data, the goal in high 
resolution photon spectrum analysis should be the suppression of noise- 
induced unphysical behaviour in the electron flux, while maintaining 
maximum ability to recover faithfully real features. Various algorithms 
have been employed with this intent (see, e.g., Johns and Lin 1992; 
Thompson et al. 1992; Piana 1994; Piana & Brown 1998: Piana et 
al. 2003) all belonging to the wide class of regularization methods 
for linear ill-posed inverse problems, though differing in method of 
regularization. For example, Johns and Lin (1992) sum over energy 
intervals to obtain sufficiently good statistical accuracy (regularization 
by coarse energy binning). Their results suggested downturn of the 
electron spectra below 50 keV, though the results were too uncertain 
to be conclusive. Piana et al. (2003) have detected, through Tiknonov 
regularized inversion (Tikhonov 1963), a feature at E ~ 55 keV in 
the mean source electron spectrum for the July 23, 2002 solar flare, 
that has been impossible to detect through a forward-fitting algorithm 
involving power-law functions such as those used by Holman et al. 
(2003). Although it has not been possible so far to establish the origin 
of this particular feature 1 , nevertheless the F(E) form obtained by 
Piana et al. (2003) is a faithful description of the F(E) corresponding 
to the photon spectrum used. 

An observed hard X-ray spectrum /(e) is related, through a bremsstrahlung 
cross-section Q(e,E), to the mean electron flux spectrum F(E) in the 
source, through the relation (Brown 1971; Johns and Lin 1992; Brown, 
Emslie & Kontar 2003) 



1 It is possible that this particular feature has non-solar origin and is a result of 
the effects of pulse pileup - Smith et al. 2002 




(1) 
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where R is the distance to the observer, V is the emitting volume and 
n = V^ 1 J n(r)dV is the mean target density. The problem of determin- 
ing F(E) from /(e) when they are related by a Volterra-type equation 
such as (1) is an ill-posed problem in the sense of Hadamard (1923). As 
a result, every experimental problem described by an equation like (1) 
is affected by a numerical pathology termed ill- conditioning whereby, 
when an unconstrained solution procedure is followed, the presence 
of measurement noise is reflected in unphysical oscillations in the re- 
constructed solution for the source function. To obtain a meaningful 
solution F(E) one needs to avoid noise amplification (e.g. Craig & 
Brown 1986) by means of regularization methods applying physical 
constraints to the electron spectrum. The algorithm looks for an ap- 
proximate least-squares solution of the integral equation relating the 
photon and the electron spectra, but subject to inclusion of additional 
information based on physically meaningful constraints or assumptions 
on the solution. This leads to the formulation of a family of regular- 
ization methods exploiting different possible a priori information on 
the electron flux coming from solar physics. The aim of the present 
paper is to introduce this generalized constraint approach into the 
field of solar HXR spectrum analysis, together with various other new 
features (physical, mathematical and numerical) yielding a more effec- 
tive algorithm, based on Generalized Singular Value Decomposition, 
for determination of F(E). More precisely, in this paper we will discuss 
the three major problems concerning the application of a regularization 
approach to the analysis of HXR data, that is: how to effectively in- 
troduce physically meaningful constraints into the inversion procedure; 
how to tune the stability requirement avoiding artificial clustering in 
the residuals, and, finally, how to adapt the regularization techniques 
to solar data characterized by huge dynamical ranges. 

The plan of the paper is as follows. In §2 we review the mathematical 
formulation of the problem pointing out its numerical instability. In 
§3 we discuss the relation between regularization and physical con- 
straints, while in §4 we describe the regularization and the computa- 
tional method based on the Generalized Singular Value Decomposition. 
In §5 a solution is proposed, through analysis of the cumulative resid- 
uals, to the crucial problem of optimal choice of the regularization 
parameter which controls the trade off between stability and loss of 
information. In §6 we perform an analysis of the solution structure. 
Finally, in §7 we present some applications in the case of simulated 
photon spectra to demonstrate the improvements achievable and the 
summary is presented in §8. 
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In a subsequent paper (Kontar et al, 2004 (hearafter Paper II)) we 
will apply this approach to real spectra observed with RHESSI (Lin et 
al, 2002). 

2. Mathematical formulation 

Equation (1) is a Volterra integral equation of the first kind and can 
be expressed in terms of a linear integral operator A : X — ► Y, 

{AF){e) = ih? ™ V l°°F( E )Q^ E ) dE , ( 2 ) 

where X and Y are two appropriate functional (Hilbert) spaces. For 
physical forms of the bremsstrahlung cross-section Q, A is a compact 
linear operator so that (Bertero et al, 1985) every discretization of 
equation (2) is characterized by (significant) numerical instabilities. 
Let us, then, consider some convenient discretized form of (1), viz. the 
linear system 

AF = g , (3) 

where 

A ij = £^Q((e i+1 + e i )/2,(E j+1 + E j )/2)5E j , i = l,...,N j = 1, . . . , M(> N)(4) 

F = (F(E 1 ),...,F(E M )), g = (s(ei),...,0(ejv)), with M > N, and 
the SEj and 5ei are appropriate weights. The values <?(ej) correspond 
to a set of discrete photon counts in energy bands — > ej + <5ej, while 
the F(Ej) are the corresponding values of the mean electron flux in 
energy bands Ej — > Ej + SEj. We use a matrix A corresponding to 
the bremsstrahlung cross-section due to Haug (1997) with the Elwert 
(1939) Coulomb correction applied 2 . 

It is important to recognize that, since electrons of all energies E > e 
contribute to the photon emission at energy e, in general the hard X-ray 
spectrum over a finite range [ei,ejv] of photon energies contains infor- 
mation on the electron spectrum over a much wider range, in particular 
within the range ejy < E < Em above the uppermost photon energy. 
For example, if F{E) has an upper energy cutoff at E = E max > ejy, 
then this cutoff will affect the observed photon spectrum below ejv, 
since the photon spectrum must tend to zero at e — ► E mSLX . Hence, by 

2 Note that the Elwert correction is not applicable for the high-frequency limit, 
e ~ E. However, since photons of energy e are produced, in general, by a wide range 
of electron energies E, the neglect of a more sophisticated form of the cross-section 
in this very narrow energy range does not significantly affect our results. 
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extending our array of E values to a sufficiently high value, the solution 
of (3) can potentially reveal evidence of an upper energy cutoff (see §7 
below) . 

Since M > N, problem (3) is under determined, i.e., there is no 
unique solution of the linear system. The same holds true if we consider 
the least-squares problem 

|| AF - g|| 2 = min (5) 
where || • || is the canonical Euclidean norm defined by 

||/|| = (/ f 2 (x)dx) (6) 

The solutions of (5) are commonly known as pseudosolutions (Bertero 
et al, 1985). Obviously, additional constraints need to be applied to 
obtain a unique solution. 



3. Physical constraints and Regularization 

As a physical quantity F must satisfy various physical conditions 
such as F > and any known constraints such as properties which are 
to be conserved or to be minimised /maximised. Many such properties 
can be expressed, for a suitable closed operator L, in the form 

||LF|| < const (7) 

leading to the need to solve least square problem (5) subject to addi- 
tional constraint (7). This constrained minimum problem can be solved 
using the Lagrange multiplier method, namely 

£(F) = ||AF-g|| 2 + A||LF|| 2 = min (8) 

where A is a Lagrange multiplier. This approach to regularize the 
problem is known as Tikhonov regularization (Tikhonov, 1963). 

Three possible different choices for L are discussed in the following. 

As a first example, we observe that the density-weighted target- 
averaged energy-integrated flux of electrons is given by the Euclidean 
norm ||F||. Physically, this quantity must be fixed by the total flux 
(the function of the total number of electrons or energy) , a requirement 
that can be met by a suitable choice of the second (constraint) term 
in equation (8), the simplest choice being ||F|| or L = 1, the identity 
matrix so that the second term of equation (8) is just the Euclidean 
norm of the solution F. Problem (8) with this constraint operator is also 
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termed zero order regularization and was used by Piana et al. (2003) 
for solar data. It defines the non-parametric F(E) of target-averaged 
total electron flux consistent with the data for a given parameter A to 
be chosen via some appropriate additional condition. 

The source averaged electron flux F(E) physically results from a 
combination of the injected electron flux Fq(Eq) spectrum and the 
physics of particle transport in the radiating source. Under a broad, but 
not comprehensive range of conditions, these flux spectra are related 
by (Brown & McKinnon, 1985) 

where dE/dN is the rate of energy loss per unit column density . For 
example, for collisional energy losses in a cold target, dE/dN ~ — 1/E 
and therefore (Brown and Emslie, 1988) 



Fo(E ) ~ -± 



F(E) 
E 



(10) 

E=E 



Equation (10) shows that, if the purpose of finding a solution F(E) 
is to subsequently use that solution to infer the injected electron flux 
spectrum Fq(Eq), then the mean electron flux should be differentiable, 
a requirement which can be incorporated in (8) by adopting for L the 
differentiation operator L ~ D 1 which is termed first order regulariza- 
tion. 

Physically, if the electron acceleration can be described determin- 
istically, e.g. the injection function resulted from acceleration can be 
presented similar to an integral (9) with systematic acceleration instead 
of deceleration, then the injected spectrum should be a differential func- 
tion. For example, acceleration by an electric field leads to differentiable 
spectrum of accelerated (injected) electrons. If we believe the injection 
function is differentiable, then F(E) should not only be differentiable, 
but have a differentiable first derivative. This corresponds to the re- 
quirement of a bounded second order derivative, hence to second order 
regularization L ~ D 2 . 

It is informative to consider the effect of applying regularization 
methods of very high order, since they can obscure physical features 
in F(E) that do not comply with the imposed smoothness constraint. 
The fc-order derivative D^F cannot be defined over less than k points, 
which reduces the dimension of the solution space to N — k for N 
data points. This makes higher order solutions more restrictive and 
potentially less precise though in practice the high resolution data 
now available have such large N ~ 300 that this is not an issue for 
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any reasonable order of regularization. However in the case of forward- 
fitting procedures, the natural iV ~ 300 dimensional space of possible 
solutions is forcefully squeezed into the space of only a few dimensions 
( 5 in the case of power-law + isothermal components - Holman et al, 
2003) which corresponds to very high order regularization and hence is 
very restrictive. 



4. Regularized solution and Generalized Singular Value 

Decomposition 

Provided that the null spaces of the matrices A and L intersect trivially 
(i.e. AF = and LF = have no common solutions other than 
F = 0), the formal solution of the minimum problem can be shown 
to be (Hansen, 1992) 

Fa = (A* A + AL*L) _1 A* . (11) 

where A* is adjoint of an operator A. From the point of view of applica- 
tions, this formula is not helpful since truncation errors imply a notable 
loss of information in forming the cross-product matrix A*A and, 
furthermore, the computational effort required is significant. Computa- 
tional heaviness, together with the presence of local minima, affect also 
the use of quadratic programming for convex functional minimization, 
which is a typical strategy for computing the solution of (8). A more 
effective approach is to use Generalized Singular Value Decomposition 
(GSVD) algorithm. Following van Loan (1976), we consider an M x N 
matrix A and a P x N matrix L (M > N > P). Then for each pair of 
real matrices (A,L) 

A G R MxN , L G R PxN . (12) 

it can be shown that there exists a set of singular values a^, satis- 
fying the relation (<J^) 2 + (cr^) 2 = 1, and singular vectors Ufc,Vfc,Wfc, 
where the first two sets are orthogonal and the third one is a set of 
linearly independent vectors satisfying the simultaneous equations 



/diag(^), 0\ 
A = U Ijv-p W -1 , L = V(diag(cr L ) 0)W' 1 . (13) 

\0 0/ 

Here the M x M matrix U is formed from the M column vectors 
Ufe, k = 1, . . . , M, with similar definitions for the P x P matrix V and 



kontar_etall.tex; 2/02/2008; 13:31; p. 7 



8 



Kontar et al. 



the N x N matrix W. The generalized singular values are defined as 
the ratios = cr^/aj^. 

The solution to this generalized minimization problem (8) can be 
shown to be (Hansen 1992) 

Fa = E ** + E (8 • (14) 

k=i\ a k + A °k ) fc=P+ i 

In the particular case when L = 1 (zero-order regularization), P = N 
and Equation (13) shows that u/% = Ufc,v^ = v&. Hence the second 
term in the solution (14) vanishes identically and the first term reduces 
to the one for the zero-order regularized solution. 



5. Choice of the Regularization Parameter 

It has to be recalled that the choice of the regularization parameter A, 
and indeed of L, has to be made independently of the equation and of 
the data, using prior knowledge or prejudice, since there is no unique 
solution to the equation (5) itself. As mentioned, the integral properties 
of the electron flux (7), if they were known, would unambiguously deter- 
mine the Lagrange multiplier A in our problem (8). Unfortunately, we 
do not know a priori the total flux of X-ray producing electrons or other 
integral quantity. Therefore it is advantageous to use knowledge of the 
errors in the recovered solution to choose the regularization parameter. 
Several criteria for determining the optimal A in Equation (14) have 
been introduced. A general property of them is that the optimal A tends 
to zero when the noise level tends to zero. For example, according to 
the discrepancy principle (Tikhonov et al, 1995), the best value of A is 
given by the solution of the equation 

||AF A -g|| 2 = |l%H 2 , (15) 

where Sg is some measure of the noise affecting the data (essentially, 
typically the canonical norm of the error vector). The discrepancy 
principle is typically rather robust in that it yields stable values but 
empirical tests show that the parameter it provides is often too large, 
and the corresponding regularized solution oversmoothed. 

Here we propose a procedure for the choice of the regularization 
parameter based on the analysis of the residuals r& = ((AF)*.— gk)/^9k- 
Then the deviation weighted by the error 

II(AFa - g)(5g) _1 || 2 ~ 1 (16) 
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accounts more accurately for point-to-point error variation than (15). 
Indeed, A denned by Equation (16) has accounted for detailed structure 
of errors, while the discrepancy principle uses only total error. 

Ideally, the normalized residuals should be consistent with statis- 
tical deviations in the data, and should therefore form a gaussian 
distribution with zero mean; the cumulative normalized residual 



should be mostly within il/y^J. A too-small set of values for the cumu- 
lative residuals indicates insufficient smoothing, whereas a set of values 
that consistently exceed ±l/y/J (especially if the sign of the residuals 
cluster) indicates too great a regularizing parameter. Therefore, we 
start with the value of A given by (16) and then reduce A until the 
average residual in the photon spectrum over the energy range from 
€\ to €j as a function of j is mostly within the ±<r limits expected if 
the residuals were purely statistical, drawn from a normal distribution. 
This technique is similar to requiring x 2 ~ 1 in hypothesis testing 
situations. 



A key issue in achieving the most meaningful construction of a regu- 
larized solution lies in analysis of the quantities 



i.e., the absolute value of the coefficients which multiply the singular 
vectors v^. in the solution (14). To obtain a meaningful solution repre- 
sented by Equation (14) the singular values should decrease faster 
on average than the coefficients (g, u^) (the Picard condition; Groetsch 
1984). Figure 1 shows a typical behaviour of the coefficients c&. 

The first singular vector is always either positive or negative definite 
(so that ciwi is always positive), while the w& always have oscillatory 
behavior for k > 1. Therefore, for sufficiently small regularization pa- 
rameters A, the regularized solution may exhibit negative values when 
the coefficients c& are not decreasing fast enough, so that relatively large 
values of the regularization parameter A are necessary to guarantee a 
positive definite solution. 

Inasmuch as typical solar flare hard X-ray spectra are sufficiently 
steep, it is advantageous to transform the fundamental problem (1) 




k=l 



(17) 



6. Analysis of the solution structure 




(18) 
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(or, equivalently, [8]) to a form in which the Ck are decreasing faster 
with k, so that smaller values of A, which preserve more fidelity in the 
recovered solution, can be used. This will also avoid errors connected 
with finite precision arithmetics in machine calculations. 

Two strategies can be followed to overcome this difficulty. In the 
first one, instead of considering equation (3) we consider the new linear 
system 

ASF = Sg (19) 

with 

<5F = F-F*; 5g = g-AF„ (20) 

where F* is an adopted form (often, but not necessarily, a closed para- 
metric expression) for F. Now Sg represents the deviations from the 
reference spectrum. Hence, if this reference spectrum is chosen appro- 
priately (say from a forward- fit solution - e.g., Holman et al. 2003), 
then these deviations SF vary around zero, and the function <5g will 
be significantly flatter than g, so that the behavior of the Ck becomes 
more monotonic. 

In addition to facilitating the calculation of a smooth (but not un- 
necessarily oversmoothed) solution of the minimization problem, the 
quantity <5F is interesting in its own right. It represents the devia- 
tion of the actual electron spectrum F(E) from the assumed reference 
spectrum F*(E) and hence is an "adjustment" to this assumed (e.g., 
forward-fitted) form. This "initial guess" approach is the one adopted 
in the present paper. 

A second possible strategy to constrain the behaviour of the c& is to 
consider the change of variables 

g(€i) -+ e? g(ei); F{Ej) -+ Ep^); A^J^A^, (21) 

with p, q positive real numbers. Then the basic form of the solution (14) 
is unaltered, but the forms of the matrix and its associated singular sys- 
tem are altered which leads to modified behavior of the coefficients c& . If 
a scaling (21) can be found that makes the decreasing functions of k, 
then this transformed solution will have more desirable properties. We 
have found through experimentation that a judicious choice of scaling 
is p = q = (7— l)/2, where 7 is the best-fit power-law spectral index to 
the array of g values. For all values of 7 (from 2 to « 20), this scaling 
drives the coefficients c\. toward a rapidly decreasing form. Note that 
while re-scaling with p = 7 produces a much flatter data function g, 
such a steeper scaling concomitantly leads (equation [21]) to the matrix 
A becoming less diagonal, thereby increasing the ill-conditioning of 
the system and hence the rate of decay of the with k. The choice 
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p = (7 — l)/2 hence lies at the "middle ground" between steepness of 
the input data vector and ill-conditioning of the transformation matrix; 
similar arguments apply to the choice of q. 

It should also be noted that "rescaling" is equivalent to constructing 
the ratio (rather than the difference) of F{E) relative to a reference 
F*(E) form. 



7. Application of the Algorithm to Simulated Data 

In this section we explore the application of the above techniques to 
simulated data, in order to demonstrate some of the effects of various 
features. 

7.1. Non-square Matrix 

To demonstrate the benefits of using a non-square A extending to 
higher electron energies than the highest observed photon energy, we 
assumed an actual F(E) of the form F(E) ~ E~ s , 10 < E < E mSLX , 
with 5 = 2 and £ max = (300,400,500) keV. We then used Equa- 
tion (1) and the electron-proton bremsstrahlung cross-section from 
Haug (1997) to calculate the photon spectrum /(e), in the range 10 < 
e < 200 keV, from each electron spectrum. (It should be noted that 
the upper limit to the "observed" photon spectrum is in all cases less 
than the upper-energy cutoff in the electron spectrum.) Figure 2 shows 
the calculated photon spectra and the corresponding local spectral 
index 7 = — dlog J(e)/c/ log e. At low energies, the spectrum is well 
represented by a power-law form with 7 = 5 + 1 = 3, but even at 
energies as low as ~ 50 keV, the deviation from a power-law behavior 
induced by the upper-energy cutoff in F(E) is already evident. For the 
case -E max = 300 keV , by e = 100 keV the spectrum has steepened 
from its low-energy form (~ e~ 3 ) sufficiently that the spectral index 
has increased by as much as 0.5. 

Such a deviation in local hard X-ray spectral index 7 is clear ev- 
idence for a significant deviation from the power-law behavior of the 
generating F(E) spectrum at higher energies, although the exact nature 
of this high-energy deviation is not immediately obvious from the form 
of the photon spectrum (or even from a 7(e) plot). We therefore now 
explore the ability of our technique to uncover the actual nature of this 
deviation (namely, in this case the upper energy cutoff at E = E mSLyi ). 
We performed a regularized inversion of the photon spectra of Fig- 
ure 2 by using solution (14) under the following conditions: zero order 
regularization; the simulated photon "data" used covered the range 
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10 < e < 200 keV and the electron upper energy limits -E U pper used in 
the inversion were 400, 500, and 600 keV. 

In all cases the recovered electron spectra quite faithfully reproduce 
the actual high-energy cutoffs E max = 300 keV (Figure 3). The higher 
energy cut-off is seen as substantial steepening in the reconstructed 
spectra. However, we should note somewhat obvious limitations: if the 
true spectra has some variations above the range we are given, the 
reconstructed solution does not display this features due to the lack of 
information given (Figure 4). 

7.2. Higher order regularization with GSVD 

To illustrate the benefits of GSVD and higher order regularization 
we consider the reconstruction of an electron spectrum which is the 
superposition of a power-law plus an oscillatory trigonometric function. 
In this case, the use of a first order penalty term governed by L ~ D 1 
is more effective than the zero-order regularization algorithm, since 
the bound on the first derivative of the regularized solution assures a 
sufficiently correct behaviour (in terms of residuals) of the solution. In 
Figure 4 (upper panel) the simulated photon spectrum is obtained by 
inserting the theoretical electron spectrum into equation (3) while Fig- 
ure 4 (lower panel) shows the reconstruction obtained by using formula 
(14) when L ~ D 1 and A is chosen by using the cumulative residuals 
analysis criterion (in this figure we also superimposed the theoretical 
electron spectrum in order to point out the reliability of our approach). 
In Figure 5(upper panel) and Figure 5(lower panel) the behaviour of the 
normalized and cumulative residuals shows that the fitting performance 
of the regularized solution is extremely accurate in the range below 150 
keV. 

The importance of higher order regularization using GSVD also 
becomes explicit when one wants to derive the injected electron dis- 
tribution Fq(Eq). Figure (6) shows both the mean source and injected 
electron spectra obtained using different orders of regularization, com- 
pared with their true forms. As a true form of Fq(Eq) we took a 
power-law with a bump simulated by exponent as can be seen in Fig- 
ure (6). The simulated photon spectrum is obtained by inserting the 
theoretical electron spectrum with the high energy cut-off at 300 keV 
into equation (3) and 5% noise has been added. The reconstructed 
spectra and input spectra are shown in the Figure (6). Note, that in 
the reconstruction we used data only up to 150 keV. 

The second order regularized solution shows the closest solution 
for Fq{Eq) in equation (6). On the other hand first and zero order 
regularizations show results with smaller x 2 in recovering the spatially 
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integrated spectra (6). Indeed, second and first order regularization 
preserves information on small scale (closer reproduction of a hump 
in Fq(Eq) around 30 keV), while zero order regularization is better in 
global properties of the solution. The first and zero order regularization 
show lack of sufficient smoothness that is displayed in oscillations in 
Fq(Eq) (6). The main deviations from the true solutions are observed 
above 200 keV, where we have only approximate solution and near low 
energy cut-off due to boundary effects. 

8. Conclusions 

In the paper, we have summarized the essential mathematics associ- 
ated with application of an Generalized Singular Value Decomposition 
(GSVD) technique to the solution of Volterra integral equations arising 
in solar X-ray spectroscopy, and in particular to the inference of mean 
source electron spectra F(E) and injected (accelerated) electron spec- 
trum Fq(Eq) from observations of solar flare hard X-ray spectra /(e). 
Judicious use of this methodology can recover forms of F(E) that are 
not only relatively free from the effects of data noise amplification, 
but which also recover features that are not realizable using more 
traditional (e.g., forward-fitting) methods. Further, they can reveal 
approximate behaviour in the electron spectrum well above the range 
of photon energies observed. 

In the companion paper (Paper II) we will make use of these tech- 
niques in the analysis of high-resolution solar flare spectra observed by 
RHESSI. 
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Figure 1. Variation of the coefficients Ck in Equation (18) as a function of the vector 
number k for the simulated data set with 8 — 2. 
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Figure 2. Upper panel: Simulated mean source electron spectra F{E), which has the 
form of a power-law (index 8 — 2), with three different upper energy cutoffs -E ma x- 
Middle panel: Corresponding photon spectra 7(e) for each of the F(E) forms. Lower 
panel: Local spectral index 7 = d log 1(e) /d log e for each photon spectrum. Note the 
evidence for the high-energy cutoff in the hard X-ray spectrum and its local spectral 
index at much lower energies than -E ma x- The vertical lines at e = 200 keV represent 
the upper energy limit of the spectral data used for subsequent analysis. 
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Figure 3. Mean source electron spectra F(E) recovered from the photon spectrum 
-E max = 300 keV) (dot-dashed spectrum in Figure 2), using the GSVD technique 
and zero order regularization with photon "data" in the range 10 < e < 150 keV and 
electron energy ranges 10 keV < E < E uppcl , where i?u PP er = 400 keV (dot-dashed 
curve), 500 keV (dashed curve) and 600 keV (dotted curve). The reference spectrum 
(Equation [20]) was of the form F,(E) ~ E 1 * . Note that the high-energy behavior 
of F(E) (in particular the high-energy cutoff at -E max = 300 keV) is quite faithfully 
reproduced. 
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Figure 4- Reconstruction from simulated data: upper panel: the simulated photon 
spectrum; bottom panel: the reconstruction obtained by using first order regular- 
ization. The dash line shows the input electron spectrum and solid lines presents 30 
realizations of the solution with the data randomly perturbed within error bars. 
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Figure 5. Reconstruction from simulated data: upper panel: normalized residuals; 
bottom panel: normalized cumulative residuals 
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Figure 6. Reconstruction from simulated data. Recovered mean flux F(E) (upper 
panel) and recovered injected flux Fo(Eo) (lower panel) given by equation (10). The 
dotted lines shows the true solution, while three various orders of regularization are 
shown: second order regularization (solid line), first order (dash line), zero (dash 
dot). 
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